Stochastic wave function method for non-Markovian quantum master equations 
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A generalization of the stochastic wave function method to quantum master equations which are 
not in Lindblad form is developed. The proposed stochastic unravelling is based on a description 
of the reduced system in a doubled Hilbert space and it is shown, that this method is capable of 
simulating quantum master equations with negative transition rates. Non-Markovian effects in the 
reduced systems dynamics can be treated within this approach by employing the time-convolutionless 
projection operator technique. This ansatz yields a systematic perturbative expansion of the reduced 
systems dynamics in the coupling strength. Several examples such as the damped Jaynes Cummings 
model and the spontaneous decay of a two-level system into a photonic band gap are discussed. The 
power as well as the limitations of the method are demonstrated. 



I. INTRODUCTION 

The theory of open quantum systems is a fundamental 
approach to the understanding of dissipation on a micro- 
scopic and macroscopic level in many fields of physics, 
such as quantum optics, and soHd state physics. Be- 
sides the conventional density matrix formaUsm, there 
has been an increasing interest over the last years in the 
stochastic wave function method [H-p|, where the state 
of the open system is described by an ensemble of pure 
states instead of a reduced density matrix. This method 
permits a description of the dynamics of an individual 
quantum system subject to a continuous measurement 
|R,p| and hence provides additional information about 
the state of the system compared to the description by 
a reduced density matrix. Moreover, the stochastic wave 
functions method has been shown to be an effective nu- 
merical tool for the solution of density matrix equations 
with many degrees of freedom since a reduced density 
matrix has 7V^ degrees of freedom, whereas a stochastic 
wave vector only has N components |§,|3 . 

Since the stochastic wave function method originated 
in the context of quantum optics most publications on 
this subject consider the weak coupling regime and the 
Born-Markov approximation seems to be inevitably con- 
nected with this approach. On the other hand, the 
stochastic wave functions method could also provide a 
useful numerical tool in other fields of physics, where 
the Born-Markov approximation is not justified. This 
has been shown for example by Imamoglu [0, who ex- 
tended the stochastic wave function method to the strong 
coupling regime by considering an enlarged system which 
contains a few ficticious bath- modes, that are weakly cou- 
pled to a Markovian environment. Although this concept 
is quite general, the drawback of this method is obvi- 
ous: even if the state space of the system is small, the 
numerical treatment of the enlarged system can become 
very expensive if more than a few ficticious modes are 
needed to approximate the reduced systems dynamics. 
Another approach to a non-Markovian stochastic wave 
function method developed by Jack, CoUett and Walls 



p2[ is based on a continuous measurement interpreta- 
tion of the stochastic unraveling. In this approach the 
stochastic equation of motion for the reduced state vec- 
tor involves a multiple time integration over the system's 
history conditioned on the measurement record over a 
finite time interval. Furthermore, it has been claimed re- 
cently by Diosi et al. |J13,H that it is in principle possi- 
ble to construct an exact stochastic Schrodinger equation 
which describes the non-Markovian time evolution of an 
open quantum system. 

In this article we present an extension of the stochastic 
wave function method beyond the Born-Markov approxi- 
mation which is based on the time-convolutionless projec- 
tion operator technique |lq;M . This ansatz yields a sys- 
tematic perturbative expansion scheme for the stochastic 
dynamics of the reduced system which is valid in an in- 
termediate coupling regime where non-Markovian effects 
are important, but a perturbative expansion is still justi- 
fied. The major advantage of our method is that it does 
not rely on an enlarged phase space and that it uses a 
stochastic evolution equation which is local in time. 

The paper is organized as follows. In Sec. || we 
briefly review three different approaches to equations of 
motion for the reduced density matrix: the derivation 
of the Markovian quantum master equation (Sec. p A| ), 
the Nakajima-Zwanzig projection operator technique 
(Sec. p^IB| ), which yields a generalized master equation, 
and the time-convolutionless projection operator tech- 
nique (Sec. II C), leading to a quantum master equation 
which is local in time. Sec. Ill deals with the stoc hastic 
unravelling of quantum master equations: in Sec. [II A 
we review the stochastic unravelling of quantum master 
equations which are in Lindblad form, such as th e Mar ko- 
vian quantum master equation, whereas in Sec. Ill B we 
present an unravelling of arbitrary linear density ma- 
trix equations which are local in time, such as the time- 
convolutionless quantum master equation. These algo- 
rithms are then applied to the damped Jaynes-Cummings 
model and the spontaneous decay of a two-level system 
into a photonic band gap in Sec. IV. Sec. V contains our 
summary. 



II. DERIVATIONS OF QUANTUM MASTER 
EQUATIONS 

We shall begin with a description of the models we 
want to examine and state the basic assumptions under- 
lying the following sections. Throughout this article, we 
consider a quantum mechanical system S which is cou- 
pled to a reservoir R. The combined system is supposed 
to be closed and its Hamiltonian is of the form 



H = Hn + aHv 



(1) 



where Hq is the Hamiltonian of the system and the reser- 
voir and Hi the interaction Hamiltonian. The time evo- 
lution of the combined system's density matrix in the 
interaction picture W{t) is determined by the Liouville- 
von Neumann equation 



d_ 
dt 



W{t) = -ia[Hi{t),W{t)] = aL{t)W{t), 



(2) 



where the interaction Hamiltonian in the interaction pic- 
ture is defined as Hi{t) = exp{iHot)Hicxp{—iHot). The 
initial state of the combined system is supposed to fac- 
torize, 



M^(0) = p(0) ® PR, 



(3) 



where pR is some stationary state of the reservoir, i. e., 
the system and the reservoir are initially uncorrelated. 
For technical simplicity, we further assume that odd mo- 
ments of Hiit) with respect to pR, vanish, i. e., 



TrR{pRffi(ii)---i7i(i2fc+i)} = 0, 



(4) 



although this assumption is not essential for the methods 
we want to use in this article (see Ref. [p^). 

Since Eq. (0) is in general a system of (infinitely) many 
differential equations, exact solutions are only known in 
rare cases. Moreover, even if an exact solution can be 
found, one is usually not interested in the dynamics of 
the environment but wants to calculate the time evolu- 
tion of system observables. Therefore, we seek an ap- 
proximate equation of motion for the reduced density 
matrix p(i) = TrR{VF(t)} of the open system. In this 
section we will describe three different approaches to that 
goal: the Born-Markov approximation, the Nakajima- 
Zwanzig projection operator technique and the time- 
convolutionless projection operator technique. 



Liouville-von Neumann equation (H) twice, differentiat- 
ing with respect to t and taking the trace over the reser- 
voir. This yields the exact equation of motion 



pW = 



-a- / dsTrR {[Hi{t), [Hi{s), Wis)]]} , 
/o 



(5) 



which still contains the density matrix W{t) of the com- 
posed system. 

The first approximation we make is the Born approxi- 
mation which consists in approximating the density ma- 
trix of the composed system by a product of the form 



VF(t) « p(i) ® PR, 



(6) 



where p{t) refers to the variables of the reduced system 
and PR denotes a stationary state of the environment. 
Such an approximation is justified if the coupling be- 
tween the system and the environment is weak. Insert- 
ing Eq. (||) into Eq. (^) we obtain the closed integro- 
differential equation for the reduced density matrix 

p(i) = -a2 / dsTrR{[iJi(i),[i?i(s),p(s)®pR]]}. (7) 
Jo 

This equation is further simplified by making the Markov 
approximation: replacing p{s) by p(i) yields a closed dif- 
ferential equation of motion for the reduced density ma- 
trix which contains only p(i) , namely 



m = -' 



dsTr 



•R{[i/i(i),[iIi(,s),p(t)®PR]]}. (8) 



The Markov approximation is based on the assumption 
that the correlation time of the reservoir tr is small com- 
pared to the time scale rs on which p{t) changes. The 
final form of the quantum master equation is obtained 
by extending the upper limit of the integral to infinity, 
which is valid for times i ^ tr since the integrand is 
negligible for s ^ tr. 

Within this derivation of the quantum master equa- 
tion, the Markov approximation appears as an additional 
approximation besides the Born approximation, and one 
is tempted to believe, that the generalized master equa- 
tion (^ is more accurate than t he in aste r equ ation (||). 
However, as we will see in Sees. II B and II C , both ap- 
proximations are only valid to second order in the cou- 
pling strength a and are hence equally accurate (see also 
|19| , p0[ ). We will also demonstrate this by means of a 
specific example in Sec. [VB, 



A. The Markovian quantum master equation 



B. Nakajima-Zwanzig projection operator technique 



In this section we sketch an intuitive derivation of the 
Markovian quantum master equation based on the Born- 
Markov approximation (see, e. g., ||l^,|l^). The start- 
ing point is an exact equation of motion for the reduced 
density matrix which can be obtained by integrating the 



The Nakajima-Zwanzig projection operator technique 
pl|-p3[ is based on a partition of the state of a system 
into a relevant and an irrelevant part by defining an ad- 
equate projection operator V which projects the state 
on the relevant part and a projector Q = 1 — V which 



projects on the irrelevant part. For our system reservoir 
model we define the projector V in the usual way as 



VW(t) = TrR {W(t)\ PR EE p(t) ® PR, 



(9) 



where pR is a stationary state of the reservoir. The equa- 
tion of motion for the two components VW{t) and QW{t) 
can be obtained directly form the Liouville-von Neumann 
equation (0): 



d_ 

di 
d_ 

dt 



VWit) = aPLit)PW{t) + aPL{t) QW{t) (10) 

QW{t) = aQL{t)rW{t) + aQL{t)QW{t). (11) 



Taking into account the initial condition Eq. (||) the for- 
mal solution of Eq. (|l^) reads 

QWit) = a [ dsg{t,s)QL{s)VW{s) (12) 

Jo 

where Q{t, s) is defined as 

g{t, s) = T^ exp (a f ds' QL{s')\ . (13) 

The symbol T^ indicates the chronological time ordering. 
Substituting the expression for QW{t) into the equation 
of motion of the relevant part of the state (|l^) we obtain 
the generalized master equation for VW{t) 

VWit)^aVL{t)VW{t)+ f dsk{t,s)VW{s), (14) 
o'l Jo 



with the memory kernel 

K{t, s) = a^VL{t)g{t, s)QL{s). 



(15) 



It is important to note, that the generalized master equa- 
tion ( |l4| ) is exact and that, hence, the explicit computa- 
tion of the memory kernel K{t, s) is, in general, as com- 
plicated as the explicit solution of the Liouville-von Neu- 
mann equation (|2|). However, Eq. (|lj) serves as a start- 
ing point for systematic approximations. For example, 
a perturbative expansion of the memory kernel K{t, s) 
to second order in the coupling strength a leads to the 
generalized quantum master equation in the Born ap- 
proximation (M). On the other hand, although the com- 
putation of the memory kernel is essentially facilitated 
by using a perturbative expansion, the final form of the 
equation of motion is still an integro-differential equa- 
tion, the integration of which can be rather difficult. We 
can overcome this by using the time-convolutionless pro- 
jection operator technique, which will be described in the 
following Section. 



C. Time-convolutionless projection operator 
technique 

The basic idea of the time-convolutionless projection 
operator technique llq,M is to replace W{s) in the for- 
mal solution of the irrelevant part (|lj) by 



W{s)=G{t,s){r + Q)W{t), 



(16) 



where the backward propagator G{t, s) of the composite 
system is defined as 



G(i, s) = T^ exp (-a f ds' L(s') j , 



(17) 



and r_^ indicates the anti-chronological time ordering. 
Solving Eq. (|l|) for QW{t), wc find 



with 



QW{t) = [1 - s(t)]"' j:{t)pw{t), 



I](i) =a f dsg{t,s)QL{s)VG{t,s) 
Jo 



(18) 



(19) 



which can be substituted in Eq. ( |10| ) to obtain the exact, 
time-convolutionless equation of motion for the relevant 
part of the system 



d_ 
dt 



VW{t) ^ K{t)VW{t) 



aVL{t)[l-Y.{t)]~^VW(t). 



(20) 



The crucial point of this construction is the existence of 
the generator K{t) which relies on the existence of the 
operator [1 — S(i)] . Since S(0) = and S(t) is con- 
tinuous, this operator exists for all t if and only if it can 
be expanded in a geometric series 



[l-S(t)]-i = ^S(t)^ 



(21) 



n=0 



This condition is always satisfied for short times or in 
the weak coupling regime, but can be violated in the 
strong coupling regime, as will be demonstrated explic- 
itly in Sec. [VB. Therefore, we define the intermediate 



coupling regime as the range of coupling parameters a, 
where non-Mar kovian effects are significant, but the gen- 
erator K{t) exists for all t. 

Using Eq. ( pl| ) we can also write the generator of the 
time-convolutionless master equation as 



K{t) = ^ aVL{t)Y.{ty 



(22) 



This form is the starting point for a perturbative expan- 
sion of K{t) in powers of the coupling strength a. To 
fourth order one obtains, for example. 



K{t) = a^K2{t) + a^Ki{t) + 0{a^) 



(23) 



where 



K2it) = / dtiVL{t)L{ti)V, 
Jo 



and 



Kiit) = / dti dt2 I dt 



2 / «t3 

^0 Jo 



(24) 



(25) 



VL{t)L{ti)L{t2)L{h)V - VL{t)L{ti)VL{t2)L{h)V 
VL{t)L{t2)VL{h)L{h)V - VL{t)L{^)VL{ti)L{t2)V 



The higher-order terms can be obtained in a way sim- 
ilar to van Kampen's cumulant expansion MJ24|. All 
terms containing odd orders of the coupling strength van- 
ish in this expansion, since by definition of V and L{t) 
we have V L{ti) ■ ■ ■ L{t2k+i)V = (sec Eq. (|)). It is 
important to note that the general structure of the time- 
convolutionless equation of motion (|2^) of the reduced 
density matrix is not changed by the perturbative expan- 
sion, i. e., the approximative equation of motion is also 
linear in p{t) and local in time, unlike the perturbative 
expansion of the generalized master equation (Oh . 



1. The equation of motion to second order 

When we substitute the expressions for the generator 
L{t) and the projection operator P, Eqs. (||) and (||), re- 
spectively, into the second order contribution to K{t) we 
immediately obtain the time-dependent quantum mas- 
ter equation (0) within the Born-Markov approximation 
(without extending the upper limit of the time integra- 
tion to infinity). Thus, Eq. (1^) as well as Eq. (0) are 
correct to the same order in the coupling and the major 
approximation in the heuristic derivation of the quantum 
master equation in Sec. HA is not the Markov, but the 
Born approximation. This seems to be somewhat coun- 
terintuitive, since after making the Born approximation 
the equation of motion of the reduced density matrix is 
still a complicated integro-differential equation, whereas 
the Markov approximation considerably simplifies the 
calculations. Nevertheless, it does in general not improve 
the accuracy of a calculation to make only the Born- 
approximation and to omit the Markov-approximation. 



2. The equation of motion to fourth order 

We now compute the explicit expression Kj^it) for 
the fourth order contribution to the time-convolutionless 
equation of motion. To this end, we decompose the in- 
teraction Hamiltonian into a sum of products in the form 



H, 



E^^- 



(26) 



TrR {pKQio {t)Qi, {ti)Qi, [t2)Q^, (is)} = 

TrR {pnQ^, it)Q,, ih)} TrR {prQ,, {t2)Q,, (is)} 
+ TrR {pKQ^o{^)Q^2 ih)} TrR {p^Q,, {ti)Q,, (is)} 
+ TrR {PKQ^MQ^. (is)} TrR {prQ,, (ii)Q,, (ia)} , (27) 

and we introduce the short-hand notation 

6,i,--- denotes i^.o(i),F,,(ti),-- • 

(12), ■ • ■ denotes TvK{pKQ^^{tl)Q^At2)} , ■ • ■ 

and sum over repeated indices ik- In this notation we 
find for example 

rL{t)---L{t3)rw{t) (28) 

= [6,[-..,[3,p]-..]](0-..3) 

= E [P^oit),[■■■AF^^{h),p{t)]■■■]] 

xTrR{pRg,„(i)---g,3(t3)}, 

and inserting the expression (Eq) into Eq. ( pq ) , we obtain 

K4{t)VW{t) ^ f dti [' dt2 [' dt3 (29) 

Jo Jo Jo 

x{(02)(13) [0, [i,2] 3p] - (02) (31) [6, [1,2] p3] 
-(20)(13) [6,3p [1,2]] + (20)(31) [6,p3 [1,2]] 
+ (03) (12) ([0, [3, 2] pi] + [0, [12, 3] p]) 
+ (30)(21)([0,ip[3,2]] + [6,p[2i,3]]) 
-(03) (21) [6, [1,3] p2\ - (30) (12) [6,2p[i,3]] }. 

Note that this expression contains commutators between 
various system operators, which can immensely simplify 
the explicit evaluation of K/^{t), if certain commutation 
relations are specified, such as bosonic commutation re- 
lations for a harmonic oscillator, or the commutation re- 
lations for the pseudospin operators (see Sec. IV A). 



III. STOCHASTIC UNRAVELLING OF 
QUANTUM MASTER EQUATIONS 

A. Quantum master equations in Lindblad form 

In Rcf. p^] Lindblad has shown that the equation of 
motion of a reduced density matrix has to be of the form 



d_ 
di 



Pit) 



^s + IE^'W^I^-z^W 



(30) 



We further assume, that the state pR is not only station- 
ary, but also Gaussian, i. e., 



E^'W ['l^l^^Pi*) - ^p{t)LlL, + L,p{t)Ll^ 



if the dynamics of the reduced system is assumed to con- 
serve positivity and to represent a quantum dynamical 



semi-group. Here, Hg is the Hamiltonian of the system, 
the time-dependent coefficients Si{t) describe an energy 
shift induced by the coupfing to the environment, namely 
the Lamb and Stark shifts, and the positive rates 7i(i) 
model the dissipative coupling to the i—th decay channel. 
In this case, the state of the open system can alter- 
natively be described by a stochastic wave function iplt) 
[y-|6| , the covariance matrix of which equals the reduced 
density matrix, i. e.. 



p{t)^ / DibDrm{^\p[^,t], 



(31) 



where P[il),t\ is the probability density functional of find- 
ing the state of the open system in the Hilbert space 
volume element D^pDip* at the time t |26|j27| ]. 

The time evolution of the stochastic wave function 
is governed by a stochastic differential equation, which 
might either be diffusive Mh or of the piecewise deter- 
ministic jump type [0-Q. The latter takes the form |0] 

d^it) = -^G{iJ,t)dt + J2(^^^^^^i,{t)^ dN^it), 

(32) 

where the dNi{t) are the differentials of indepen- 
dent Poisson process Ni{t) with mean {dNi{t)) — 
'yi{t)\\Liip{t)\\'^dt. The drift generator takes the form 



G{i^,t) = H{t)i; + -Y^S.mlL.^I; 



2^-^ 



T.^m{Ll 



L, 



|i»^f W 



(33) 



For the differential of the Poisson process dNi (t) the Ito 
rule dN,{t)dNj{t) = S^jdNi{t) holds, that is, dN,{t) can 
either be or 1. If dNiit) = 0, then the system evolves 
continuously according to the nonlinear Schrodinger-type 
equation 

i-^{t) = G{i,,t), (34) 

whereas, if dNi(t) — 1 for some i, then the system un- 
dergoes an instantaneous, discontinuous transition of the 
form 



^{t) 



L,^{t) 



(35) 



Note that the generator G'(V', t) of the continuous time 
evolution is non-Hermitian and hence the propagator of 
ipit) is non-unitary. However, due to the nonlinearity of 
the generator, the norm of i/;(i) is preserved in time. 

Using the Ito calculus for the differentials dNi(t) it is 
easy to check that the equation of motion of the covari- 
ance matrix of ipjt) equals the usual Markovian quantum 
master equation ( |3CJ ) in Lindblad form. Thus, expecta- 
tion values of system observables can either be calculated 
by means of the reduced density matrix or as averages 
over different realizations of the stochastic process ip{t) 
and both descriptions yield the same results. 



B. General quantum master equations 



The most general type of a quantum master equation 
which results from the time-convolutionless projection 
operator technique - or from a perturbative approxima- 
tion - is linear in p{t) and local in time (see Sec. pC) but 
needs not to be in the Lindblad form, as we will show in 



an example below (see Sec. IVC). However, these equa- 
tions can always be written in the form 



-Pit) = Ait)pit) 



pit)B\t) + J2ait)pit)DUt), 



(36) 



with some time-dependent linear operators A{t), B{t), 
Ci{t), and Di(t). In order to find an unraveling of this 
equation of motion we follow a strategy, which has al- 
ready been successfully applied to the calculation of mul- 
titime correlation functions |^,^ . We describe the state 
of the open system by a pair of stochastic wave functions 



9it) = 



( m 



Ht) 



(37) 



Formally, 6{t) can be regarded as an element of the dou- 
bled Hilbert space n = H ® H. If P[9,t] denotes the 
probability density functional of the process in the dou- 
bled Hilbert space 7i, we may define the reduced density 
matrix as 



p(i)- DdD9*\<l>){^\P[e,t] 



(38) 



The time evolution of the state vector 9{t) is then gov- 
erned by the stochastic differential equation 



d0(t) ^-iG{0,t)dt 



(39) 



■E 



II^WII 



WMmt) 



-Mt)e{t)-9{t)]dN,{t) 



where dNi (t) is the differential of a Poisson process with 
mean 



(dN-(t)) = M(M^d, 

and the functional G{9,t) is defined as 

l^||J.(t)^(i)| 



(40) 



Gie,t)=^[F{t) + -J2 



2^ \\0{tW 
with the time-dependent operators 



e{t), (41) 



, _ ( A{t) 
^^^i" y B{t) 



Mt) 









A(i) 



(42) 



Again, this type of stochastic evolution equation de- 
scribes a piecewise deterministic jump process, where the 
deterministic pieces are solutions of the differential equa- 
tion 



i^^e{t)^G{e,t), 



(43) 



and the jumps induce transitions of the form 

Oi^)^m)Ljm^J^(?^i]. (44) 



11^.^(01 



\\mt)\ 



Note, that the structure of the stochastic differential 
equation in the doubled Hilbert space ( p9| ) is very simi- 
lar to the structure of the stochastic differential equation 
(p2|). In fact, the unraveling of general quantum master 
equations presented in this section contains as a special 
case the u nraveling of Lindblad-type equations shown in 
Sec. |[IIA|: If we set 



A{t) = B{t) - -iHs -IY. [^^(^) + *^'^(^)] 4^fe (45) 



and 



a{t) = D,{t) ^ ./^L,, 



(46) 



the equation of motion (^6|) reduces to the Lindblad equa- 
tion (5fl) and both unravelings are identical. 



IV. EXAMPLE: THE SPONTANEOUS DECAY OF 
A TWO-LEVEL SYSTEM 



In this section we consider as an example of the general 
theory the exactly solvable model of a two-level system 
spontaneously decaying into the vacuum within the ro- 
tating wave approximation. The Hamiltonian of the total 
system is given by 

Ho = wscr+cr" + J2 ^kblbk, (47) 

k 

iJi = CT+ B + cr" ® fit ^ith B^Y^ cjkhk, (48) 



where ujs denotes the transition frequency of the two- 
level system, the index k labels the different field modes 
with frequency Wfe, annihilation operator bu and coupling 
constant gk, and ct^ denote the pseudospin operators. 



A. Exact and approximated equations of motions 

The exact solution and equation of motion for this 
model can be obtained in the following way: Define the 
states [M 



V'o = |0)s ® |0)r 

V'i = |i)s®|o)r 

^k = |0)s® |fc)R 



(49) 



where |0)s and |l)s indicate the ground and excited state 
of the system, respectively, the state |0)r denotes the vac- 
uum state of the reservoir, and |fc)R = &|,|0)r denotes the 
state with one photon in mode k. Since the interaction 
Hamiltonian conserves the total number of particles, the 
flow of the Schrodinger equation generated by Hi is con- 
flncd to the subspace spanned by these vectors. Hence, 
we may expand the state of the total system at any time 
as 



0(t) = coV'o -f ci(i)f/'i +^Cfc(i)V'fc, 



(50) 



with some probability amplitudes cq, ci(t), and Ck{t). 
The time evolution of these probability amplitudes is de- 
termined by a complicated system of ordinary differential 
equations, which can be solved in some simple cases by 
introducing the so-called pseudomodes [gO|. With these 
probability amplitudes, the reduced density matrix takes 
the form 



Pit) 



|ci(t)|2 Ci(t)cS 

ct(t)co |coP + EJcfe(t)p 



(51) 



Differentiating this expression with respect to time we 
get the following exact equation of motion. 



lpit) = -'-S{t)[a+a-,p{t)] 



(52) 



-lit) 



1 



a+a p{t) p(f)o-+(T" + a~p(t)a^ 



where the time-dependent energy shift S(t) and decay 
rate 7(t) are deflned as 



S{t) = -23 



cijt) 
ci{t) 



lit) 



'2» 



ciit) 
ciit) 



(53) 



Note, that if the decay rate 7(i) is positive for all t, then 
this equation of motion is in the Lindblad form ( po[ ) . 

The equation of motion within the Born approxima- 
tion can be expressed in terms of the reservoir correlation 
function. To this end, we define the real functions <i>(t) 
and *(t) as 



$(t) + i^{t) = 2TrR {B{t)B^pR} e^"=^* 
2 / dLoJ{u;)e'^'^^-'^^\ 



(54) 



where B(t) = e'xjp{iHot)B e:cp{—iHot), and we have per- 
formed the continuum limit. J(w) is the spectral density. 



The equation of motion in the Born approximation (m) 
then reads 






(55) 



-cr+cr p{s) + -p{s)a+a -a p{s)a'^ 



Performing the Markov approximation and extending the 
upper hmit of the time integral to infinity, we obtain the 
usual time-independent quantum master equation 



u t 



(56) 



where the Markovian Lamb shift Sm and the Markovian 
decay rate 7m are defined as 



5, 



M 



ds^{s), 7m 



ds<^{s) 



(57) 



The time-convolutionless expansion of the equation of 
motion according to Sec. II C leads to a quantum master 
equation which has the same structure as the exact equa- 
tion of motion, but the time-dependent energy shift S{t) 
and decay rate 7(i) are approximated by the quantities 

5(4) (t) = / dh^it -ti) + l f dh f dt2 f ' dt3 

Jo ^ Jo Jo Jo 

X ^{t - t2)$(ii - h) + $(< - ^2)^(^1 - ta) 



*(t - t3)^ti - t2) + $(t - t3)*(ti - t2) 



and 



7(4) (t) = / dii$(t-ti 
Jo 



1 



dti 



dU 



JO Jo 

*(i - t2)^{h - h) - $(f - t2)$(ti - h) 

*(t - i3)*(ti - t2) - Ht - hMh - h) 



(58) 



dU 



(59) 



It is important to note that the explicit expressions for 
S"'"*) (t) and 7''*' (t) only involve ordinary integrations over 
the reservoir correlation functions, which can be done an- 
alytically in simple cases or numerically. 



B. The damped Jaynes-Cummings model on 
resonance 

The damped Jaynes Cummings model describes the 
coupling of a two-level atom to a single cavity mode which 
in turn is coupled to a reservoir consisting of harmonic 
oscillators in the vacuum state. If we restrict ourselves to 
the case of a single excitation in the atom-cavity system. 



the cavity mode can be eliminated in favor of an effective 
spectral density of the form 



Jiu) 



1 



7oA2 



27r (ws - ^)2 + A2 ' 



(60) 



where cug is the transition frequency of the two-level sys- 
tem. The parameter A defines the spectral width of the 
coupling, which is connected to the reservoir correlation 
time tr by the relation tr ~ A^^ and the time scale 
Ts on which the state of the system changes is given by 
Ts = lo^- "^^"^ exact probability amplitude ci(t) (see 
Eg pO|)) is readily obtained by using the method of poles 
||30| , since J{uj) has simple poles at tu — utQ zL iX. One 
gets 

ci(t) = ci(0)e-^*/2 ("coshl + ^sinh^) , (61) 



where d = \/X? — 270 A, which yields the time-dependent 
population of the excited state 



Pii{t) 



dt 



pii(0)e-^* (cosh 2 



A , dt 
— smh — - 
d 2 



(62) 



Using Eq. (p3) we therefore obtain a vanishing Lamb 
shift, S{t) = 0, and the time-dependent decay rate 



lit) 



270A sinh(di/2) 



(63) 



dcosh{dt/2) + Asinh(dt/2) ' 

In Fig. n (a) we illustrate this time-dependent decay rate 
7(t) ('exact') together with the Markovian decay rate 
7m = 7o ('Markov') for rg = 5rR. Note that for short 
times, i. e., for times of the order of tr, the exact decay 
rate grows linearly in t, which leads to the quantum- 
mechanically correct short -time behavior of the transi- 
tion probability. In the long-time limit the decay rate 
saturates at a value larger than the Markovian decay 
rate, which represents corrections to the golden rule. The 
population of the excited state is depicted in Fig. ^ (b) : 
for short times, the exact population decreases quadrati- 
cally and is larger than the Markovian population, which 
is simply given by /3ii(0) exp (—70^), whereas in the long- 
time limit the exact population is slightly less than the 
Markovian population. 

Next, we want to determine the solution of the gener- 
alized quantum master equation in the Born approxima- 
tion. To this end, we insert the spectral density of the 
couphng strength (|^ into Eq. (|4|) to obtain *(t) = 
and 



(64) 



$(i)==7oAexp(-At). 



The solution of the generalized master equation (pq) can 
be found in the following way. We differentiate Eq. (|55| ) 
with respect to t and obtain 



m = -\p{t) 



70 A 



(65) 



-(j+(j~p{t) p(t)cr+(T~ + a~p{t)a+ 



Due to the exponential memory kernel, this equation of 
motion is an ordinary differential equation which is local 
in time, and contains only p{t), p{t) and p{t). Solving 
this system of differential equations for p(i), we obtain 
the time evolution of the population of the upper level 



Pii(t) = pii(o)e-^*/2 L^^^ Y^i ^'""^ t) 



(66) 



where d' = y/X^ — 470A. From this expression, we can 
determine the time-dependent decay rate 



m 



27oAsinh(d't/2) 



Pii(t) _ 

pii (t) d' cosh{d't/2) + A sinh(d'i/2) ' 



(67) 



the structure of which is similar to the exact decay rate 
(|63|). Note, however, the difference between the parame- 
ters d and d' which can also be seen in Fig H (a) where 
we have also plotted the decay rate 7(i) ('GME 2'): For 
short times, the decay rate j(t) is in good agreement with 
7(t), but in the long time limit, "/{t) is too large. 

Finally, the time-convolutionless decay rate can be de- 
termined from Eq. (p9), and to second and fourth order 
in the coupling we obtain 



7(^^(0 = 70 (1-e-^*) 



(68) 



and 



7W(t)=^o 1 



,-At 



To 
A 



sinh(At)-At]e-^*}, (69) 



respectively, which corresponds to a Taylor expansion 
of the exact decay rate j{t) in powers of 70, as can be 
checked by differentiating 7(t) with respect to 70. Fig. |l| 
(a) clearly shows, that j^^^'it) as well as 7^'*-'(i) approx- 
imate the exact decay rate very good for short times, 
and 7*-'*^(i) is also a good approximation in the long time 
limit. 

The time evolution of the population of the excited 
state can be obtained by integrating the rate 7*^''' (t) with 
respect to t. This yields 



pitHO=Pii(0)exp 



ds-r^-^Hs) 



(70) 



In order to compare the quality of the different approx- 
imation schemes, we show the difference between the 
approximated populations and the exact population in 
Fig. |] (c). Besides the analytical solutions of the gener- 
alized master equation (pq) and the time-convolutionless 
master equations, we have also performed a stochas- 
tic simulation of the time-convolutionless quantum mas- 
ter equations with 10^ realizations. Since the approx- 
imated rates j^^'^^t) are positive for all i, the corre- 
sponding master equations are in Lindblad form, and 
we can use th e stochastic simulation algorithm described 
in Sec. Ill A as an unravelling. Fig. |l| (c) shows, that 
the stochastic simulation is in very good agreement with 
the corresponding analytical solutions. Moreover, we 



see that the difference between the solution of the time- 
convolutionless master equation to fourth order and the 
exact master equation is small (see also Fig || (b)), 
whereas the errors of the generalized and the time- 
convolutionless master equation to second order which 
correspond to the Born approximation and the Born- 
Markov approximation (without extending the integral) , 
respectively, are larger and of the same order of magni- 
tude. In fact, the Markov approximation even leads to 
a slight improvement of the accuracy, compared to the 
Born approximation, which is surprising if we consider 
the heuristic derivation of the quantum master equation 
in Sec. |IIA| . 

As we pointed out in Sec. O, the approximation 
schemes used in this article are perturbative and hence 
rely on the assumption that the coupling is not too 
strong. But what happens, if the system approaches the 
strong coupling regime? We will investigate this question 
by means of the damped Jaynes-Cummings model on res- 
onance, where the explicit expressions of the quantities 
of interest are known. 

First, let us take a look at the exact expression for 
the population of the excited state (p2|): In the strong 
coupling regime, i. e., for 70 > A/2 or rs < 2rR, the pa- 
rameter d is purely imaginary. Defining d = —id we can 
write the exact population as 



pii(i) = pii(0)e 



-At 



dt X . dt\ ,_ , 



which is an oscillating function that has discrete zeros at 

(72) 



2 / d 

t — — \ TT ~ arctan — 

d \ A 



Hence, the rate 7(i) diverges at these points (see 
Eq. ([5^)). Obviously, 7(t) can only be an analytical func- 
tion for t G [O,io[i where to is the smallest positive zero 
oipii{t). 

On the other hand, as we have seen in Sec. [VB , 
the time-convolutionless quantum master equation corre- 
sponds basically to a Taylor expansion of j(t) in powers 
of 7o, and the radius of convergence of this series is given 
by the region of analyticity of j{t). For 70 < A/2, this is 
the whole positive real axis, but for 70 > A/2 the pertur- 
bative expansion only converges for t < to. This behavior 
can be clearly seen in Fig. n^ (d), where we have depicted 

pii{t) and Pii (t) for rs = tr/5, i. e., for a strong cou- 
pling: the perturbative expansion converges to pii(i) for 
t '^S to ~ 6.3/70, but fails to converge for t > to. 

The solution of the generalized master equation to sec- 
ond order shows a quite distinct behavior, but also fails 
in the strong coupling regime: for 70 > A/4 the popu- 
lation pii{t) starts to oscillate and even takes negative 
values, which is unphysical (see Fig. n] (d)). 

The 'failure' of the time-convolutionless master equa- 
tion ai t — to can also be understood from a more in- 
tuitive point of view. The time-convolutionless equation 



of motion ( |2C)| ) states that the time-evolution of p{t) only 
depends on the actual value of p{t) and on the gener- 
ator K{t). However, at t = to the time evolution also 
depends on p{0). This fact can be seen in Fig. g, where 
we have plotted pii{t) for three different initial condi- 
tions, namely pii{0) — 1.0, 0.5, 0.0. At t — to, the cor- 
responding density matrices coincide, regardless of the 
initial condition. However, the future time evolution for 
t > to is different for these trajectories. It is therefore 
intuitively clear that a time-convolutionless form of the 
equation of motion which is local in time ceases to exist 
for t > to. The formal reason for this fact is that at t = to 
the operator 1 — E(t) (see Sec. HC), is not invertible and 



hence the generator K{t) does not exist at this point. 



C. The damped Jaynes-Cummings model with 
detuning 

In this section we treat the damped Jaynes-Cum mings 
model with detuning, i. e., the same setup as in Sec. IV B 
but the center frequency of the cavity coq is detuned by 
an amount A = us—ujo against the atomic transition fre- 
quency. In this case the spectral density of the coupling 
strength reads 



J(c.) 



1 



7oA2 



and thus the functions $(t) and ^(t) are given by 

$(t) =7oAe"^*cos(At), 
*(t) =7oAe"^*sin(At). 



(73) 



(74) 
(75) 



With these functions, the time-dependent Lamb shift and 
decay rate to fourth order in the coupling, S''^-'(t) and 



X '{t), respectively, can be calculated using Eqs. (|5^ 
and (p9|). The integrals can be evaluated exactly and 
lead to the expressions 



A2-^A2 

7gA2A3e-^* 
"2(A2 + A2)3 



{[l-3(ir]( 



.At 



e"^* cos 



(2At)) 



l-(^) Atsin(At)+4 l + (^) Atcos(At) 



3-(^)' 



e~^* sm 



in(2At)} 



(76) 



and 



l''\t) 



loX' 



A2+A2 



[l-e-^*(cos(At)-f sin(At))] 



7o^A^e 



-xt 



2(A2 



A2)3 
,4 



{[i-3(fr]( 



-,At 



e~^* cos 



i{2At)) 



1_(A)'' Atcos(At)+4 l+(^) Atsin(At) 



3 -(f) 



e '^* sm 



in(2At)|. 



(77) 



In Fig. ^ (a) we have depicted 7^'*''(t) together with the 
exact decay rate, which can be calculated using the meth- 
ods of poles [0 for A = 8A and A = 0.870. Note, that 
the spontaneous decay rate is severely suppressed com- 
pared to the spontaneous decay on resonance. This can 
also be seen by computing the Markovian decay rate 7m 
which is given by 



7M 



7oA2 



A2 



0.01570. 



(78) 



However, this strong suppression is most effective in the 
long time limit. For short times, 7(t) oscillates with a 
large amplitude and can even take negative values, which 
leads to an increasing population. This is due to photons 
which have been emitted by the atom and reabsorbed at a 
later time. Hence, the exact quantum master equation as 
well as its time-convolutionless approximation are not in 
the Lindblad form (BOh , but conserve the positivity of the 
reduced density matrix! This is of course not a contradic- 
tion to the Lindblad theorem, since a basic assumption 
of the Lindblad theorem is that the reduced system dy- 
namics constitutes a 1-parameter dynamical semi-group. 
However, in our example this is not the case, since the 
initial preparation singles out the specific time t = and 
the domain of the operator K{t) is shrinking for increas- 
ing t. 

Since the transition rate 7'-^-' (t) also takes negative val- 
ues, we can not use th e stochastic simulation algorithm 
presented in Sec. HI A for a stochastic unravelling of the 



time-convolutionless quantum master equation, but have 
to use the simu lation algorithm in the doubled Hilbert 
space (see Sec. |IIIB ). The dynamics of the stochastic 
wave function 9{t) = (0(t), V'(t))"'", which is an element 
of the doubled Hilbert space H = H®H is governed by 
the stochastic differential equation (|39|), where the oper- 
ators F and J are given by 



F = -lj(^\t) 



a ' (7 








and 



J 



|7(4)| 



sign (7'-'*') cr 




(79) 



(80) 



The deterministic part of the time evolution is governed 
by the nonlinear Schrodinger-type equation 



|.(.)^G,M)^= F+i^ .W, (81, 



which results in a continuous drift, whereas the jumps 
induce instantaneous transitions of the form 



0(t) 



\\jem 



je{t) 



y(4)) 

|0)s 



sign ( 7^'M |0)s 



(82) 



If the rate j^^\t) is positive then this type of transition 
leads to a positive contribution to the ground state popu- 
lation poo{t): whereas a negative rate leads to a decrease 
of Poo- 

In Fig. (b), we show the results of a stochastic simula- 
tion for 10^ realizations, together with the analytical so- 
lution of the time-convolutionless quantum master equa- 
tion and the exact solution. Obviously, the agreement of 
all three curves is good and the stochastic simulation al- 
gorithm works excellently even for negative decay rates. 
In addition, we also show the solution of the Markovian 
quantum master equation which clearly underestimates 
the decay for short times and does not show oscillations. 



D. Spontaneous decay into a photonic band gap 

As our final example, we treat a simple model for the 
spontaneous decay of a two-level system in a photonic 
band gap which was introduced by Garraway Q. To 
this end, we consider a spectral density of the coupling 
strength of the form 



V. SUMMARY 



JH = 



"o 
2tt 



WiTi 



(c.-c.s)2 + (ri/2)2 



(^-C.s)2 + (r2/2)2/' ^^^^ 

where ilg describes the overall coupling strength, Fi the 
bandwidth of the 'flat' background continuum, F2 the 
width of the gap, and Wi and 14^2 the relative strength 
of the background and the gap, respectively. Again, the 
function J{u)) has a small number of poles, and hence 
the exact solution can be determined by using pseudo- 
modes [gfj. In Fig. ^ (a) we show the excited state's 
decay rate j{t) for the same parameters as in Ref. [ pi[ , 
i. e., Fi/^o = 10, F2/flo = 1, VFi = 1.1, and 1^2 = 0.1. 
For short times, 'y{t) increases linearly on a time scale of 
Fj^ and then takes a maximum, which stems from tran- 
sitions into the 'flat' background continuum. For longer 
times, i.e., i 3> Fj"^, the transitions into the background 
are suppressed, and the decay rate becomes smaller and 
smaller until it reaches its final value. Thus, the popula- 
tion of the excited state decreases rapidly for times of the 
order F^^, and slowly in the long-time limit (see Fig. 

(b)). 

The time-dependent Lamb shift S^"^' (t) and the decay 
rate ^^'^'{t) of the time-convolutionless quantum master 
equation to fourth order can be computed by inserting 
the spectral density of the coupling strength J(w) into 
Eq. (II). This yields ^(i) = and 



<^{t) = 2nl(Wi 



-rit/2 



Woe 



-r2t/2 



(84) 



which can be inserted into Eqs. ( |5§| ) and (p9|). Since 
*(i) = the Lamb shift S^^\t) vanishes; the time- 
dependent decay rate 7^*-'(i) can be computed explicitly, 
and is in good agreement with the exact decay rate for 
our choice of parameters (see Fig. || (a) and (b)). 



In this article we have presented a generalization of 
the stochastic wave function method to quantum master 
equations which are not in Lindblad form. This general- 
ization - together with the use of the time-convolutionless 
projection operator technique - makes it possible to ex- 
tend the range of potential applications of the stochas- 
tic wave functions method beyond the weak coupling 
regime, where the Born-Markov approximation is valid, 
without enlarging the system. This generalization is ca- 
pable of treating systems in the intermediate coupling 
regime, i. e., systems for which the generator of the time- 
convolutionless quantum master equation exists for all t 
and is analytic in the coupling strength a. In the exam- 
ples we investigated in this article, this range was limited 
by Ts > Tpj. The dynamics of this class of systems is gov- 
erned by an equation of motion which is local in time 
and can be approximated by a perturbative expansion. 
This perturbative expansion leads in general to a quan- 
tum master equation, which needs not to be in Lindblad 
form but can be unravelled with our method. The basic 
idea of this unravelling is the introduction of stochas- 
tic processes in a doubled Hilbert space, which has been 
already successfully used for the computation of matrix 
elements of operators in the Heisenberg picture and mul- 
titime correlation functions. 
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FIG. 1. The damped Jaynes-Cummings model on resonance. Exact solution (exact), time-convolutionless master equation 
to second (TCL 2) and fourth order (TCL 4), generalized master equation to second order (GME 2), and the Markovian 
quantum master equation (Markov): (a) Decay rate of the excited state population, (b) the population of the excited state, 
including a stochastic simulation of the time-convolutionless quantum master equation with 10^ realizations, and (c) deviation 
of the approximative solutions from the exact result, for rs = Str (moderate coupling), (d) Population of the excited state for 
rs — 0.2rH, (strong coupling). 
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FIG. 2. The damped Jaynes-Cummings model on resonance. 
Pii(O) = 1.0, 0.5, 0.0 in the strong coupling regime (rs — 0.2tr). 



Exact population for the three different initial conditions 



0.15 
0.10 



° 0.05 



to 



0.00 



-0.05 



-0.10 








exact 

TCL4 

Markov 



12 



18 



70* 




(a) -^ (b) 

FIG. 3. The damped Jaynes-Cummings model with detuning. Exact solution (exact), time-convolutionless master equation 
to fourth order (TCL 4), and the Markovian quantum master equation (Markovian): (a) Decay rate of the excited state 
population, and (b) the population of the excited state, including a stochastic simulation of the time-convolutionless quantum 
master equation with 10^ realizations, for A = 0.870 and A = 8A 
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FIG. 4. Spontaneous decay in a photonic band gap: Exact solution (exact), and time-convolutionless master equation to 
fourth order (TCL 4): (a) Decay rate of the excited state population, and (b) the population of the excited state, including 
a stochastic simulation of the time-convolutionless quantum master equation with 10^ realizations, for Wi — 1.1, 14^2 ~ 0.1, 
Ti/fio = 10, and Ta/no = 1. 
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